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■ Abstract - We demonstrate that the first arrival time in dispersive processes in self-affine fractures 
O , are governed by the same length scale characterizing the fractures as that which controls their 

■ permeability. In one-dimensional channel flow this length scale is the aperture of the bottle neck, 
i.e., the region having the smallest aperture. In two dimensions, the concept of a bottle neck is 

, generalized to that of a minimal path normal to the flow. The length scale is then the average 

aperture along this path. There is a linear relationship between the first arrival time and this 
length scale, even when there is strong overlap between the fracture surfaces creating areas with 
zero permeability. We express the first arrival time directly in terms of the permeability. 
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. Due to their role in the flow properties of tight and low permeability reservoirs such as 

' shale gas reservoirs and carbonate reservoirs, and on contaminant transport e.g. in connec- 

■ tion with waste storage, the study of transport in fractures is still a very vigorous field [IHl] ■ 
Most present theoretical efforts attempts to relate the transport properties of fractures to 
the statistics of the aperture fields through analytical models based on statistical averages, 
weak disorder perturbation expansions [5] , mean- field approximations or simplified aperture 
models |B]. We also mention the work of Zhan and Yortsos [7] where a method to deduce 

' the heterogenities of a permeability distribution from the concentration arrival time field 

was proposed. 

Due to the surface roughness, i.e., the heterogenities of the aperture field, these relations 
provide satisfactory results only over a finite range of conditions and do not permit to 
predict the behavior of a fracture with large heterogenities in aperture field. One of the 
main difficulties is to correctly take into account the increasing influence of the contact 
area as the fracture aperture is decreased fSHTT]. We analyze in this Letter the dispersion 
problem at finite Peclet number and identify the proper aperture measure for this problem, 
taking into account severe heterogeneities such as large contact zones. Our main focus is on 
the breakthrough time, i.e., the time at which the tracer appears at a given position. The 
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Fig. 1: Fluid flows in a self-affine aperture field seen from above. The curves are the streamlines 
of the flow field. The streamlines have been found using the Kirchhoff method. We have then used 
the Lattice Boltzmann method to simulate dispersion fixing the Peclet number at 10. The tracer 
concentration is shown on a grey scale where darker means higher concentration. Areas where the 
aperture is zero — i.e., the fracture surfaces are in contact — are shown as black. 

surprising result that we find is that this aperture is the same as the one controlling the 
permeability P^[T^ . 

There are now numerous experimental studies and field observations that demonstrate 
that natural fractures have a self-afiine roughness |14f[T5] — for a review, see Bonamy and 
Bouchaud [T2| . Self-affine fractures are characterized by a scaling invariance of the statistical 
properties of the surface roughness under a rescaling of the distances by a factor A in the 
average fracture plane and a rescaling of the heights. Here C, is the Hurst or roughness 
exponent which takes value close to 0.8 for rocks such as granite [20] and values close to 0.5 
for porous rocks like sandstone pTll22] . 

Wc consider in this work synthetic self-affine fracture surfaces that have been generated 
using a Fourier method [inilll]- The fracture is modeled by matching the fracture surface 
with an opposite flat surface. Since we use the Reynolds approximation, this correctly 
models flow in fractures as it is only the aperure that enters the Sow equations. We define 
the fracture aperture H(r) as the distance between the two surfaces at position r. In 
the present work, the rough surface progressively approaches the fiat surface one and the 
aperture of overlapping regions is set to zero. Hence, H[r) > where there is no overlap 
and H{r) = where there is overlap. 

The flow field is determined for a fixed pressure difference between the fracture inlet and 
outlet by solving the finite differenced Reynolds equations through LU decomposition. The 
total volume entering the fracture per width and time at the inlet is proportional to the 
pressure AP over the fracture, and is given by 



(1) 



where K is the permeability and n the viscosity of the fluid. L is the length of the fracture. 
The breakthrough time is typically obtained by summing along each streamline the 
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local time of convection [25H27] . This method does, however, not take into account diffusion 
between and along the streamlines. We have instead discretized the velocity field on a square 
lattice with nearest-neighbor and next nearest-neighbor connections. If we assume that the 
dispersion time between neighboring nodes is the inverse of the velocity component along the 
vector between them, we have then introduced diffusion into the dispersion process. This is 
so since a given tracer particle will not follow the streamlines but move between the nodes 
via projected velocity vectors. We add to this description the analysis of Stern [28] of the 
first arrival time of a diffusive process with or against a convective velocity field, making it 
possible to tune the Peclet number. In order to find the first arrival time, we use the optimal 
path algorithm of Hansen and Kertesz [53] . We have verified our algorithm by comparing it 
to a two-dimensional Lattice Boltzmann method [50] • Fig. [T] shows the tracer concentration 
in gray levels at breakthrough with Peclet number Pe — 10 based on the Lattice Boltzmann 
method. In the following, we do not discuss any dependence of our results on the Peclet 
number. The numerical experiments we report have been done at Pe = 10 as a reasonable 
value. Other values add nothing significant. 

Transport properties of fractures are often characterized by equivalent or apparent aper- 
tures — such as the hydraulic, mass balance or electrical apertures — which refers to the 
aperture of a fracture with flat and parallel walls having the same property as the original 
fracture. In practice, equivalent apertures are estimated from hydraulic and conservative 
tracer tests. The mass balance aperture bm is defined as the ratio between the fluid flux Q 
and the averaged fluid velocity u [51115^ . In practice, the average fluid velocity u equals the 
average of the velocities of all fluid particles and should be derived from the average resi- 
dence time determined from the momentents of the time distribution of the measured tracer 
breakthrough curve |31| . Here we connect the first arrival time, Tmim to the average fluid 
velocity u by the expression u — LjTmin, as proposed by Guimera and Carrera [32]. Since 
the pressure gradient is kept fixed, the mass balance aperture hm = Q/u is proportional to 
Kt 

Before considering two-dimensional fractures, i.e., fractures where the aperture is or- 
thogonal to a two-dimensional fracture plane, we consider a one- dimensional version of the 
problem, i.e., a fracture where the aperture is orthogonal to a fracture line. We introduce 
a cartesian coordinate system with the x axis along the line which now consitutes the flat 
surface. Let us set a = minx H{x). We then define 

h{x) = H{x) - a. (2) 

When a < 0, the fracture is closed and hence the permeability is zero. For positive a, it is this 
parameter that controls the permeability of the fracture in the lubrication approximation 
[l3] . The permeability is in this limit given by the expression 

K io kic^C+af ' ^-"^ 

where k and c are two parameters, c, the topothcsis which characterizes the roughness of 
the aperture field, is a length scale, and k has the units of permeability. For large a, this 
gives rise to the scaling relation 

K ^ La^ , (4) 

whereas for intermediate a, we find 

K ^ La^-^/^ . (5) 

For small a, the permeability is completely controlled by the region where h{x) = 0, and 
the continuum approach behind eq. @ breaks down. We then find that the permeability is 
given by 

K - iOa^ . (6) 
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We now calculate T,nin using the lubrication approximation and in the infinite Peclet 
number hmit where diffusion is absent. The first tracer to traverse the rough channel is the 
one which has traveled along the streamline located midway between the walls i.e. where 
the velocity is at its maximum. The time this has taken is Tmm, and it is given by 



^ dx 
u{x) ' 



(7) 



where u{x) is the maximal velocity at position x along the channel which is proportional to 
the flow rate over the local aperture in the lubrication aproximation. Wc have thus 



1 



Tmin « 7^ / {K^) +a)dx . (8) 

y JO 



Combining eqs.((T]) and (|5]), we get 

T,„i„i^=^^ {h{x)+a)dx. (9) 

This integral may be performed by using order statistics [13]: We order the function h{x) 
h[^] — h(x\^ such that < /i[^2] when < ^2- For a sclf-affine profile we have that 

~ and eq. ([5]) becomes 

where c is a constant. Hence, we have the central result for a one-dimensional channel 

T^inK = A + Ca, (11) 

where A oc and C (x L^. Hence, we see that it is the minimum aperture a which 

controls the first arrival time Tmin. This is the same aperture that controls the permeability, 
see eqs. Q - ([S]). This is a somewhat surprising result, since at the minimal apperture 
location, because of mass conservation, the flow rate is maximal. Consequently, the time in 
the bottle neck effect has a small contribution to the integral eq. ([7]). However, the bottle 
neck controls the total flow rate, and this, in turn, controls the first arrival time. 

As described in Talon et al. [12j[T3], the extrapolation of the bottle neck effect to a 
two-dimensional fracture is not straight forward since the minimal aperture point is easily 
bypassed by the flow. However, it is possible to generalize the concept to two dimensions by 
replacing the minimal aperture a by the minimal path aperture. In order to introduce this 
concept, we orient our fracture such that the one of the edges parallel to the average flow 
direction follows the x axis. The y axis follows the edge where the tracer is injected and 
the z axis is orthogonal to the average fracture plane. Hence, < a; < i and < y < W^. 
We define C{x) as a path starting at (x, y = 0) and ending at (x' , y = W) without crossing 
itself. Hence, we define the quantity 



B{x) ^ 



WL 



min / di-ey{lY 

C(^) Jc(x) 



1/3 

(12) 



This is the minimal average fracture opening over all paths starting at (x, 0) and ending 
anywhere along the opposite edge at y = W . This quantity corresponds to H{x) in the 
one- dimensional case. The minimal path aperture is defined as 

6c = minS(x), (13) 
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Fig. 2: Loglog plot of TminK — A as function of a and be for C, = 0.3 (upper) and ^ = 0.8 (lower). 
The solid curves are for r-mi-nK — A vs. be, whereas the broken curve is for Ti^i-nK — A vs. a. We 
determine A by varying it until we obtain the best possible power law. In both figures, the straight 
portions of the curves have unit slope as indicated in eq. (|19|) . The curves are based on one sample 
of size 512 x 512 for each roughness. 



corresponding to the smallest aperture a in the one-dimensional case. We finally define 

b{x) = B{x) - 6e , (14) 

in the same way we defined h{x) in eq. ([2]) in the one-dimensional case. 

The central idea in Talon et al. [13] was that the two definitions be and h{x) could replace 
a and h{x) in the one-dimensional case in the permeability integral ([3]). After ordering 
h{x) ->■ b[^], we find that b[^] - ^1^, where P = 1.5 for C = 0.8 and /3 = 1.2 for C = 0.3. The 
intermediate scaling regime (O then is replaced by 

K - WLbl-^//^ , (15) 

whereas the large and small scale regimes become respectively 

K - WLbl , (16) 

and 

K - WL°bl . (17) 
Numerical experiments based on solving the Kirchhoff equations give 

r for C = 0.8, 

^ ^ 1 VFL62 i6±o.o2 fo^. ^ ^ 0.3 , 



for the intermediate regime. The results are very close to the prediction of eq. ([T5|) . 

By following exactly the same procedure for the first arrival time, i.e., replace h(x) by 
b{x) and a by be in eq. ([5]), we find 

T^inK = A + Cbe, (19) 
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Fig. 3: Scaling of coefficients A (upper figure) and C (lower figure) defined in eq. (|f9|l . The straight 
lines are best fits and have slopes 3 for the = 0.8 data (+) and 2.2 for the = 0.3 data (x) in 
the upper figure. The theoretical values are 2.8 and 2.3 respectively. In the lower figure, the best 
fit has slope 2.0 for both the = 0.8 and = 0.3 data. The theoretical value is 2. 

where A oc WL"^^^ and C oc WL^ . Fig. [3] shows A and C as a function of L verifying these 
two scaling laws. We show in fig. [21 Tmin-K' as a function of he and of a. We see that the 
linearity of T,^{nK is verified for the entire range of he values, whereas it is only true for large 
values of a. This is where a and he begin to coincide. 

Hence, we have verified that the miminal path aperture he controls both the permeability 
K and the minimal time r^m^. This, together with cq. (jl9p constitute two main results of 
this Letter. 

Since the first arrival time and the permeability are controlled by the same aperture, 
it is possible to eliminate the aperture between them. Hence, we may express the first 
arrival time directly in terms of the permeability by combining eqs. (jl6p - (|19p. We show 
in Fig. |3]T,„in vs. K for two roughnesses, C = 0.8 and C = 0.3. We expect that for small A', 
Tmin ~ A/ K'^ ^ whcrc A is defined in cq. (jl9p . For large /<, we expect Tmin '-^ C/K'^^^. For 
intermediate K, we expect T,„i„ - A/K + C/K^-^'^ when C = 0.8 and Tmin ~ A/K + C/K^-^'^ 
for C = 0.3. C is defined in eq. ((T9)) . Straight lines with the appropriate slopes have been 
added in Fig. 21 For the intermediate region, we have used only the term proportional to C . 

We have in this Letter shown that the permeability and the first arrival time in dispersive 
processes are controlled by the same aperture length scale in self-affine fractures. We have 
also shown that the functional relation between the first arrival time and this aperture 
is linear, see eq. (jl9p . The appropriate aperture is the minimal path aperture defined in 
eqs. (HH) and It is a generalization of the concept of the narrowest constriction that 

controls both the permeability and the first arrival time in one-dimensional fracture systems. 
Whereas the scaling properties we report are specific to self-affine aperture fields, the method 
of analyis based on optimal paths is not. 

* * * 

We thank J. P. Hulin for many interesting discussions. A. H. thanks the Universite 
de Paris-Sud 11 for financial support. H. A. and L. T. thanks the PICS "The Physics 
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Fig. 4: First arrival time Tmin vs. permeability K for C, = 0.8 (upper figure) and C, — 0.3 (lower 
figure). Each curve is based on one sample of size 512 x 512. The slopes of the straight lines are 
given in the main text. 

of Geological Complex System" and the Rcseaux de Thcmatiques de Rcchcrches Avancees 
"Triangle de la Physique" for financial support. 
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